{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# 3.6 DG/HDG splitting\n",
    "\n",
    "When solving unsteady problems with an operator-splitting method it might be benefitial to consider different space discretizations for different operators.\n",
    "\n",
    "For a problem of the form\n",
    "\n",
    "$$\n",
    "  \\partial_t u + A u + C u = 0\n",
    "$$\n",
    "\n",
    "We consider the operator splitting:\n",
    "\n",
    "* 1.Step: Given $u^0$, solve $t^n \\to t^{n+1}$: $\\quad \\partial_t u + C u = 0 \\Rightarrow u^{\\frac12}$\n",
    "* 2.Step: Given $u^{\\frac12}$, solve $t^n \\to t^{n+1}$: $\\quad \\partial_t u + A u = 0 \\Rightarrow u^{1}$\n",
    "\n",
    "This splitting is only first order accurate but allows to take different time discretizations for the substeps 1 and 2. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In this example we consider the Navier-Stokes problem again (cf. [unit 3.2](../unit-3.2-navierstokes/navierstokes.ipynb)) and want to combine \n",
    "\n",
    "* an $H(div)$ -conforming Hybrid DG method (which is a very good discretization for Stokes-type problems) and\n",
    "* a standard **upwind DG** method for the discretization of the convection\n",
    "\n",
    "The weak form is: Find $(\\mathbf{u},p):[0,T] \\to (H_{0,D}^1)^d \\times L^2$, s.t.\n",
    "\\begin{align}\n",
    "\\int_{\\Omega} \\partial_t \\mathbf{u} \\cdot v + \\int_{\\Omega} \\nu \\nabla \\mathbf{u} \\nabla \\mathbf{v} + \\mathbf{u} \\cdot \\nabla \\mathbf{u} \\ \\mathbf{v} - \\int_{\\Omega} \\operatorname{div}(\\mathbf{v}) p &= \\int f \\mathbf{v}  && \\forall \\mathbf{v} \\in (H_{0,D}^1)^d, \\\\ \n",
    "- \\int_{\\Omega} \\operatorname{div}(\\mathbf{u}) q &= 0 && \\forall q \\in L^2, \\\\\n",
    "\\quad \\mathbf{u}(t=0) & = \\mathbf{u}_0\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Again, we consider the benchmark setup from http://www.featflow.de/en/benchmarks/cfdbenchmarking/flow/dfg_benchmark2_re100.html . The geometry:\n",
    "\n",
    "![](../unit-3.2-navierstokes/geometry.png)\n",
    "\n",
    "The viscosity is set to $\\nu = 10^{-3}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from netgen import gui\n",
    "from ngsolve import *\n",
    "from netgen.geom2d import SplineGeometry\n",
    "geo = SplineGeometry()\n",
    "geo.AddRectangle((0, 0), (2, 0.41), bcs = (\"wall\", \"outlet\", \"wall\", \"inlet\"))\n",
    "geo.AddCircle((0.2, 0.2), r=0.05, leftdomain = 0, rightdomain = 1, bc = \"cyl\")\n",
    "mesh = Mesh( geo.GenerateMesh(maxh = 0.08) ); Draw(mesh)\n",
    "order = 3\n",
    "mesh.Curve(order)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For the HDG formulation we use the product space with\n",
    "\n",
    "* $BDM_k$: $H(div)$ conforming FE space (degree k)\n",
    "* Vector facet space: facet functions of degree k (vector valued and only in tangential direction)\n",
    "* piecewise polynomials up to degree $k-1$ for the pressure\n",
    "\n",
    "![](resources/stokeshdghdiv.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## HDG spaces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "V1 = HDiv ( mesh, order = order, dirichlet = \"wall|cyl|inlet\" )\n",
    "V2 = FESpace ( \"vectorfacet\", mesh, order = order, \n",
    "              dirichlet = \"wall|cyl|inlet\" )\n",
    "Q = L2( mesh, order = order-1)\n",
    "V = FESpace ([V1,V2,Q])\n",
    "\n",
    "u, uhat, p = V.TrialFunction()\n",
    "v, vhat, q = V.TestFunction()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parameter and directions (normal / tangential projection):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "nu = 0.001 # viscosity\n",
    "n = specialcf.normal(mesh.dim)\n",
    "h = specialcf.mesh_size\n",
    "def tang(vec):\n",
    "    return vec - (vec*n)*n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Stokes discretization\n",
    "The bilinear form to the HDG (see [unit-2.8-DG](unit-2.8-DG/DG.ipynb) for scalar HDG) discretized Stokes problem is:\n",
    "\\begin{align*}\n",
    "K_h((\\mathbf{u}_T,\\mathbf{u}_F,p),(\\mathbf{v}_T,\\mathbf{v}_F,q) := & \\displaystyle \\sum_{T\\in\\mathcal{T}} \\left\\{ \\int_{T} \\nu {\\nabla} {\\mathbf{u}_T} \\! : \\! {\\nabla} {\\mathbf{v}_T} \\ d {x} - \\int_{\\Omega} \\operatorname{div}(\\mathbf{u}) q \\ d {x} - \\int_{\\Omega} \\operatorname{div}(\\mathbf{v}) p \\ d {x} \\right. \\\\\n",
    "& \\qquad \\left. - \\int_{\\partial T} \\nu \\frac{\\partial {\\mathbf{u}_T}}{\\partial {n} }  [ \\mathbf{v} ]_t \\ d {s}\n",
    "- \\int_{\\partial T} \\nu \\frac{\\partial {\\mathbf{v}_T}}{\\partial {n} }  [ \\mathbf{u} ]_t \\ d {s}\n",
    "+ \\int_{\\partial T} \\nu \\frac{\\alpha k^2}{h} [\\mathbf{u}]_t [\\mathbf{v}]_t \\ d {s}  \\right\\}\n",
    "\\end{align*}\n",
    "\n",
    "where $[\\cdot]_t$ is the tangential projection of the jump $(\\cdot)_T - (\\cdot)_F$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "alpha = 4  # SIP parameter\n",
    "dS = dx(element_boundary=True)\n",
    "a = BilinearForm ( V, symmetric=True)\n",
    "a += nu * InnerProduct ( Grad(u), Grad(v)) * dx\n",
    "a += nu * InnerProduct ( Grad(u)*n,  tang(vhat-v) ) * dS \n",
    "a += nu * InnerProduct ( Grad(v)*n,  tang(uhat-u) ) * dS\n",
    "a += nu * alpha*order*order/h * InnerProduct(tang(vhat-v),  tang(uhat-u))*dS\n",
    "a += (-div(u)*q -div(v)*p) * dx\n",
    "a.Assemble()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The mass matrix is simply\n",
    "\n",
    "\\begin{align*}\n",
    "M_h((\\mathbf{u}_T,\\mathbf{u}_F,p),(\\mathbf{v}_T,\\mathbf{v}_F,q) := \\int_{\\Omega} \\mathbf{u}_T \\cdot \\mathbf{v}_T dx\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "m = BilinearForm(V , symmetric=True)\n",
    "m += SymbolicBFI( u * v )\n",
    "m.Assemble()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Initial conditions\n",
    "As initial condition we solve the Stokes problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "gfu = GridFunction(V)\n",
    "\n",
    "U0 = 1.5\n",
    "uin = CoefficientFunction( (U0*4*y*(0.41-y)/(0.41*0.41),0) )\n",
    "gfu.components[0].Set(uin, definedon=mesh.Boundaries(\"inlet\"))\n",
    "\n",
    "invstokes = a.mat.Inverse(V.FreeDofs(), inverse=\"umfpack\")\n",
    "gfu.vec.data += invstokes @ -a.mat * gfu.vec\n",
    "\n",
    "Draw( gfu.components[0], mesh, \"velocity\" )\n",
    "Draw( gfu.components[2], mesh, \"pressure\" )\n",
    "Draw( Norm(gfu.components[0]), mesh, \"absvel(hdiv)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Application of convection\n",
    "\n",
    "In the operator splitted approach we want to apply only operator applications for the convection part.\n",
    "Further, we want to do this in a usual DG setting. \n",
    "As a model problem we use the following procedure:\n",
    "\n",
    "* Given $(\\mathbf{u},p)$ in HDG space: project into $\\hat{\\mathbf{u}}$ in usual DG space\n",
    "* Solve $\\partial_t \\hat{\\mathbf{u}} = C \\hat{\\mathbf{u}}$ by explicit scheme (involves convection evaluations and mass matrix operations only)\n",
    "* Solve Unsteady Stokes step with r.h.s. from convection sub-problem. To this end evaluate mixed mass matrix $\\int \\hat{\\mathbf{u}} \\cdot \\mathbf{v}$ to obtain a functional on the HDG space"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For the projection steps we use mixed mass matrices:\n",
    "\n",
    "* $M_m$ : $HDG \\times DG \\to \\mathbb{R}$\n",
    "* $M_m^T$ : $DG \\times HDG \\to \\mathbb{R}$\n",
    "* $M_{DG}$ : $DG \\times DG \\to \\mathbb{R}$ (block diagonal)\n",
    "\n",
    "To set up mixed mass matrices we use a bilinear form with two different FESpaces. \n",
    "\n",
    "We do not assemble the matrices as we will only need the matrix-vector applications of $M_m$, $M_m^T$ and $M_{DG}^{-1}$.\n",
    "\n",
    "There is a more elegant version (`ConvertL2Operator`, similar to `TraceOperator`) that we discuss [below](#convertl2)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### mixed mass matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "VL2 = VectorL2(mesh, order=order, piola=True)\n",
    "uL2, vL2 = VL2.TnT()\n",
    "bfmixed = BilinearForm ( V, VL2, nonassemble=True )\n",
    "bfmixed += vL2*u*dx\n",
    "bfmixedT = BilinearForm ( VL2, V, nonassemble=True)\n",
    "bfmixedT += uL2*v*dx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### convection operator\n",
    "* convection operation with standard Upwinding\n",
    "* No set up of the matrix (only interested in operator applications)\n",
    "* for the advection velocity we use the $H(div)$-conforming velocity (which is pointwise divergence free)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "vel = gfu.components[0]\n",
    "convL2 = BilinearForm(VL2, nonassemble=True )\n",
    "convL2 += (-InnerProduct(Grad(vL2) * vel, uL2)) * dx\n",
    "un = InnerProduct(vel,n)\n",
    "upwindL2 = IfPos(un, un*uL2, un*uL2.Other(bnd=uin))\n",
    "\n",
    "dskel_inner  = dx(skeleton=True)\n",
    "dskel_bound  = ds(skeleton=True)\n",
    "\n",
    "convL2 += InnerProduct (upwindL2, vL2-vL2.Other()) * dskel_inner\n",
    "convL2 += InnerProduct (upwindL2, vL2) * dskel_bound"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Solution of convection steps as an operator (`BaseMatrix`)\n",
    "We now define the solution of the convection problem for an initial data $u$ in the HDG space. The return value (\"res\") is $M_m \\hat{u}$ where $\\hat{u}$ is the solution of several explicit Euler steps of the convection problem\n",
    "$$\n",
    "  \\partial_t \\hat{u} = - M_{DG}^{-1} C \\hat{u}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "class PropagateConvection(BaseMatrix):\n",
    "    def __init__(self,tau,steps):\n",
    "        super(PropagateConvection, self).__init__()\n",
    "        self.tau = tau; self.steps = steps\n",
    "        self.h = V.ndof; self.w = V.ndof # operator domain and range\n",
    "        self.mL2 = VL2.Mass(Id(mesh.dim)); self.invmL2 = self.mL2.Inverse()\n",
    "        self.vecL2 = bfmixed.mat.CreateColVector() # temp vector\n",
    "    def Mult(self, x, y):\n",
    "        self.vecL2.data = self.invmL2 @ bfmixed.mat * x # <- project from Hdiv to L2\n",
    "        for i in range(self.steps):\n",
    "            self.vecL2.data -= self.tau/self.steps * self.invmL2 @ convL2.mat * self.vecL2\n",
    "        y.data = bfmixedT.mat * self.vecL2\n",
    "    def CreateColVector(self):\n",
    "        return CreateVVector(self.h)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The operator splitted method (Yanenko splitting)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Setup of $M^*$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "t = 0; tend = 0\n",
    "tau = 0.01; substeps = 10\n",
    "\n",
    "mstar = m.mat.CreateMatrix()\n",
    "mstar.AsVector().data = m.mat.AsVector() + tau * a.mat.AsVector()\n",
    "inv = mstar.Inverse(V.FreeDofs(), inverse=\"umfpack\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The time loop:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "tend += 1\n",
    "res = gfu.vec.CreateVector()\n",
    "convpropagator = PropagateConvection(tau,substeps)\n",
    "while t < tend:\n",
    "    gfu.vec.data += inv @ (convpropagator - mstar) * gfu.vec\n",
    "    t += tau\n",
    "    print (\"\\r  t =\", t, end=\"\")\n",
    "    Redraw(blocking=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## `ConvertL2Operator` <a class=\"anchor\" id=\"convertl2\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Instead of doing the conversion between the spaces \"by hand\", we can use the convenience operator `ConvertL2Operator` that realizes\n",
    "$ M_{DG}^{-1} \\cdot M_m $:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "HDG_to_L2 = V1.ConvertL2Operator(VL2) @ Embedding(V.ndof, V.Range(0)).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "The transpose of $ M_{DG}^{-1} \\cdot M_m $ is $ M_m^T \\cdot M_{DG}^{-1}$\n",
    "* `V1.ConvertL2Operator(VL2)`: V1 $\\to$ VL2\n",
    "* `V1.ConvertL2Operator(VL2).T`: VL2$'$ $\\to$ V1$'$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Overwrite the application and use the `ConvertL2Operator` instead ot the mixed mass matrices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "def NewMult(self, x, y):\n",
    "    self.vecL2.data = HDG_to_L2 * x # <- project from Hdiv to L2\n",
    "    for i in range(self.steps):\n",
    "        self.vecL2.data -= tau/self.steps*self.invmL2@convL2.mat*self.vecL2\n",
    "    y.data = HDG_to_L2.T @ self.mL2 * self.vecL2    \n",
    "    \n",
    "PropagateConvection.Mult = NewMult"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The time loop:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "tend += 1\n",
    "convpropagator = PropagateConvection(tau,substeps)\n",
    "while t < tend:\n",
    "    gfu.vec.data += inv @ (convpropagator - mstar) * gfu.vec\n",
    "    t += tau\n",
    "    print (\"\\r  t =\", t, end=\"\")\n",
    "    Redraw(blocking=True)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
